Recursive quantum algorithm to find the lowest eigenstate of a general Hamiltonian 
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We propose a recursive quantum algorithm to find the lowest eigenstate of a general Hamiltonian. 
It yields the lowest eigenstate directly from an arbitrary chosen initial state without diagonalizing 
the Hamiltonian matrix. Remarkably, the effectiveness of the algorithm does not depend on any 
specific Hamiltonian structure as far as the initial state is not orthogonal to the lowest eigenstate. 
Further, our algorithm does not get trapped in the local-minima of the Hamiltonian. The number 
of recursions required for high accuracy ~ 1 - e (e « 1) are bounded by the order of 
where D is the difference between the two lowest eigenvalues. 

PACS numbers: 03.67.Ac, 02.60.Dc 



Dealing with complex systems whose dynamics are 
typically described by high (possibly a few hundred or 
thousand) dimensional Hamiltonian is essential in vari- 
ous field of modern science In particular, one of the 
frequently faced problems is how to find the eigenstate 
of the Hamiltonian associated with the lowest eigenvalue 
@~B|]- Solving this problem is of theoretical and practi- 
cal importance in many problems, such as studies of spin 
glasses @, neural netwroks [f|, and protein structures 
[3, H[ . It is also related to the global optimization such 
as planning @ and scheduling [l(J. Most problems of 
this class belong to intractable problems (so-called NP- 
complete) due to the cost of computational power as high 
as that required for the diagonalization of whole Hamil- 
tonian matrix (even when only the lowest few eigenstates 
are needed); the number of eigenstates s oug ht are inde- 
pendent on the difficulty of the problem |Tl| . 

Several methods have been proposed to find the lowest 
eigenstate 0, [lll - flEj , but their large overheads in dealing 
with a high dimensional Hamiltonian is still a detrimental 
factor in practical applications. Various efforts have been 
devoted so far to improve these algorithms [Tfj| - |20j , or to 
develop novel algorithms with a reduced calculation time 
[3, [2lH23| . However, in most methods, the system Hamil- 
tonian is typically assumed to have a specific structure: 
diagonal elements are dominant in the Hamiltonian ma- 
trix (diagonally dominant) [U, [25| , or only few elements 
distributed irregularly throughout the matrix are nonzero 
(sparse) [TTL [261] . Although these restrictions on Hamilto- 
nian would be valid and useful in several problems e.g. in 
ground-state calculations for molecules, solids, and sur- 
faces [27j , or in the case when the computational bases to 
characterize the sparseness of the Hamiltonian are known 
for each particle, an alternative effective approach may be 
essential in regard of more general Hamiltonian problems 
in complex systems composed of many particles (28|. 

Moreover, in some previous iterative methods, a trial 
wave vector typically generated at random is often signifi- 
cantly employed (e.g. in density-matrix renormalization- 
group method [27]) so that an inappropriately chosen 



trial wave vector possibly cause the algorithm to be 
trapped in local minima of the Hamiltonian 2!) ii{It . It 
can thus result in a waste of computational resources. 

In this paper, we propose a recursive quantum algo- 
rithm to find the lowest eigenstate of a general Hamil- 
tonian. In our algorithm, we first prepare an arbitrary 
chosen initial state in the same dimensional Hilbert space 
and construct a sufficiently large, but finite, number of 
transformation operators in a recursive way. Then, by 
applying the transformation operators to the initial state, 
we can obtain a state close to the lowest eigenstate un- 
less the initial state would not be in a state orthogonal 
to the lowest eigenstate. We note that our algorithm is 
generally applicable to the case when the Hamiltonian 
is neither diagonally dominant nor sparse. Further, it 
does not get trapped in local minima. We show that the 
number of recursions required for high accuracy ~ 1 — e 
(e <C 1) are bounded by the order of 0((De) _1 ), where 
D is the difference between the two lowest eigenvalues. 

Let us start with a brief overview of our algorithm. 
We consider a Hamiltonian H in iV-dimensional Hilbert- 
space, having the eigenvalues A& scaled in order such that 
< Aq < Ai < • • • < Ajv_i < 1, where no perfect degen- 
eracy between Ao and Ai is assumed, i.e., Ai — Ao ^ 0. 
In this circumstance, our goal is to find the lowest eigen- 
state | Ao) - To begin, we prepare an arbitrary pure ini- 
tial state IV'o) in -/V-dimensional Hilbert-space. We then 
construct a sufficiently large (but finite) number n c of 
transformations, • • ■ , T nal based on self reference, 

i.e., a transformation of the current step is made by that 
of earlier steps. In this sense, our algorithm is recur- 
sive. Applying the transformations to \ipo) in turn, we 
can obtain a desired output state \tpn a )- 



T nc --- T 2 Tx |V>o) = \lpn e ) 



(1) 



which is close to | Ao) - Therefore, constructions of Tj's 
(i = 1, 2, • • • , n c ) are at the heart of our algorithm. 

From now on, we shall describe details of the pro- 
cedure of our algorithm. We first choose an arbitrary 



2 



state |<jf>o) in A^-dimensional Hilbert-space, represented 
as \(j)o) — Ylj=o a j \ w j) w ith computational bases |iOj-) 
for j = 0, 1, • • • , N — 1. Here it is assumed that the coef- 
ficients aj and \wj) are completely known to us. We also 
assume that \(f>o) is not orthogonal to the lowest eigen- 
state (4>o\Xq) ^ 0. We then use a single-qubit ancillary 
system initiated as the diagonal state (|0) + \ l))/y/2 to 
improve the stability of our algorithm by removing the 
fluctuation of the probability to obtain the lowest eigen- 
state. Thus, the composite initial state is 




(2) 



Arbitrarily state 



If we expand the state |<^o) in terms of the eigenstates 
|Afc) of H, it is rewritten by 



|V>o) - 



V2 



N-l 



E a i ( A fcK) l A fe) 



N-l 



^^(|0,A fe } + |l,A fe )) ; 



fe=0 



(3) 



where |0,A fc ) = |0) ® \X k ), |l,A fe ) = |1) ® |A fe ), and 
7o.fe = SjLo 1 a i i^k\wj). Note that 70^ cannot be com- 
puted because the eigenstates |Afc) are unknown to us. 

In order to construct the transformations, we define a 
2iV-dimensional unitary operator with the Hamiltonian 
H such that 

U = (|0) (0| ® + i (|1) (1| ® i(/s) f ) , (4) 

where A(k) = e lK,H indicates a unitary transformation 
in the iV-dimensional Hilbert-space, and K is a scaling 
factor concerning the time of Hamiltonian action. The 
scaling factor k is chosen in (0, 7r/3]. This condition is 
important for successful application of the algorithm, as 
described in later. 

Then, we introduce another useful operator, so-called 
Householder reflection (HR), which is widely used in nu- 
merical linear algebra [2J, [32[ . The standard form of HR 
is written as, 



R a = 1 - 2 la) (a\ 



(5) 



where |a) is a state and 1 is identity. This can be un- 
derstood as a reflection about a state vector \a) with its 
relative phases |33f | . 

We now define two element transformations used in 
our recursive algorithm: For i = 1, 2, • • • , n c , 



Ri = TiRi-iT- , 



(6) 



where Rq — 1 — 2 |-0o) (ipo\- Thus, starting from Rq 
[I — 2 |^o) (V'ol with the (known) initial state \ipo), we 
can construct all transformations without any reliance 
on unknown or artificially introduced variables. Such a 
bottom-up construction is enabled when a process can be 
separated into sub-processes whose results can be used in 
a larger process. This methodology is known to save the 
computational power in calculations [34l [35| . 

The transformations are applied on \tpo) as in Eq. (jTJ) . 
Then the initial state \ipo) evolves progressively: 



IV>1 



T2, 



1-02 



r 3 , 



(7) 



where the states \ipi) is the output state of ith recursion 
step defined as 



\ipi) = Ti |V»— 1) = 2iTi. 



' ?1 |Vo 



(8) 



The state can be rewritten in terms of the eigen- 
states, for i — 1,2, ... , n Cl such that 



N-l 



iVi) = e t| ( n Ti,* 10, A fc > + n it, 1 1 • -V' ■ ) • ( " > 



fe=0 



Here the coefficient 7^ (or its complex conjugate 7* fc ) is 
given as a function of fcth-ordered eigenvalue A& : 



7i,fc = (4W?_i - 1) - 2Wi_ l£ 
where the factor is defined by 



-ik(1— Afe) 



(Vil^lVi 



(10) 



(11) 



Observing Eq. ||5J), Eq. (fTO)) and Eq. (|TT|). we can see that 
the coefficients 7^ are also given in a recursive manner. 
Note that they also cannot be computed due to the un- 
known eigenstates |Afe). 

Given \ip ne ), a projection is performed on the ancillary 
qubit onto 0/1 basis to remove the state of the ancillary 
qubit. The probability of the remaining state being kth- 
ordered eigenstate \X k ) is then given by 



Pn c (A/s) = |7n c ,fc| 2 l7n c -i,fe| 2 • • • l7i,fe| 2 po(Afe) 

= |7n c ,fc| 2 Pn c -l(Afe), 



(12) 



where po(\k) = |7o,fc| ■ Here, we remind that our al- 
gorithm is applicable for the case when (0o|Ao) 7^ 0, or 
equivalently, pq(Xq) 7^ 0. We note that the probability 
of picking up randomly the initial state \<po) orthogonal 
to I Ao) is extremely small, and our algorithm works well 
even with extremely small overlap between them. 

What we shall show is the proof that the probability 
Pi(Ao) of the lowest eigenstate increases with increasing 
the recursion i. Firstly, for convenience, we define |7i,fc| 2 
in Eq. (fT2|) as "quotient" such that 



Pi(X k ) (Xk\ [Tr Afipi^f J )\X k ) 
>i,k = rrr = 7 \ 1 U-j) 



(A* I (' 



Tr A /5i_i I A, 



3 



where Tta denotes partial trace over the ancillary qubit 
and pi = \ipi) Next, we write the factor Wj in 

Eq. HI]) as 



N-l 



Wj = ^ Pt(Afc) cos k (1 - A fe ) 



From Eq. (| 10[) and Eq. p^|) . we can see that 



!n c ,0 



> Qn c ,l > • • • > Qn c: 



JV-1- 



(14) 



(15) 



Then, as the scaling factor k is a constant chosen in 
(0, 7r/3], the factor Wj is bounded as 



i < Wi < 1, 



(16) 



where the lower bound is given in the limit of p;(Ao) — > 1 
with Pi(Xk) — (k ^ 0) and Ao — 0; the upper bound 
is given in the limit of pi(XN-i) — > 1 with pi(Xk) — 
(k^ N-l) and Ajv-1 ~ 1. Then, by using Eq. dTUJ) and 
Eq. (fTB"]) . we can also prove that 



Qi,o > 



> 



> Qr 



> 1. 



(17) 



From the Eq. (TT51) and Eq. (IT71) , it can be concluded that 
the probability Pi(Xq) of the lowest eigenstate is always 
increased as the recursion proceeds increasing i as 



Po(Ao) < • • • < p„ c _i(Ao) < p™ c (A ). 



(18) 



We note that our algorithm does not get trapped in the 
local-minima of the Hamiltonian, which is distinct from 
other standard search algorithms |36j |. 

Finally, we shall show that Pi(Arj) approaches ~ 1 by 
finite number n c of the recursions. This is very important 
as it is closely related to the efficiency of the algorithm in 
the use of computation power. To do this, firstly, suppose 
that we have p„ c _i(Ao) = 1 — e, with very small error e. 
From Eq. (jlll) . we can estimate W nc -i as below: 



W nc 



(1 — e) cos k (1 — Ao) 



k(1-Ai) 



cos ft + kAq sin/t + eft |Ai — Aq| shift, (19) 



where Ao and Ai are assumed to be small so that the 
terms with higher order of them are negligible. Such an 
estimation is consistent with the above analysis. Insert- 
ing Eq. (Ti"9l) mto Eq. (fTOj) . the quotient Qn c ,o at n c th step 
is given, up to the first order of e (vanishing e 2 — > 0) as 



Qn c ,o — 1 + 2eK (sin 4k + sin 2k) D, 



(20) 



where D = Ai — Ao is the difference between the two low- 
est eigenvalues. Then, by using Eq. ([15)) andp„ c _i(Ao) = 
1 — e, the increment of the probability is approximately 
given by 

PnJXo) ~ p„ c _i(A ) + 2en (sin 4k + sin2/t) D. (21) 

Therefore, Pi(Xo) continuously approaches ~ 1 for very 
small e, and equal to 1 when e = 0. From Eq. (|2ip . wc 
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FIG. 1: (a) We draw p<(Ao), p<(Ai), and p<(A 2 ) for AT = 10 3 
with respect to i. We assume a specific whose eigenvalues 
are given as Eq. (I23|l . The scaling factor is chosen to be 
7r/4. For simplicity, | Vo) is prepared with 70, k = l/\fN. As 
analyzed, po(^o) grows to ~ 1. (b) The graph of Q iy k versus i 
is also given. The data are also consistent with our analysis. 



can also prove that the recursion number n c required for 
an accuracy of ~ 1 — e is bounded as 



n c < O 



1 

7b 



(22) 



which is given in terms of D and e for very large N . 

Let us now investigate the effectiveness of our al- 
gorithm with a specific example. We consider a A^- 
dimensional Hamiltonian H with equidistant eigenvalues 
such that 



Afc = 



fc+1 



(fc = 0,l,--- ,N-1), 



(23) 



which is simple but known to be significant [37H39j . We 
apply our algorithm for N = 10 3 and plot Pi(Xo), Pi(Xi), 
and Pi(X2) in Fig. [TJa) as increasing the number of re- 
cursion i. We set the scaling factor k = 7r/4, and, for 
simplicity, |i/>o) is prepared with 70^ = 1/y/N for all k. 
It shows that pi(Xo) increases up to ~ 1, whereas Pi(Xi) 
and Pi(Xz) decrease down to zero. We also plot the quo- 
tients Qi t k (k = 0,1,2) against the number of recursion 
i in Fig. [Ub). It is observed that Q i is always larger 
than 1. On the other hand, Qi i and Qi^ become smaller 
than 1, resulting in decreases of the probabilities Pi(X\) 
and Pi(X 2 ). 

We also investigate n c against N and e as plotted in 
Fig. [21(a) on a log-log scale. We compare three differ- 
ent cases where |Ai — Arj| is given to be l/N, 1/N 2 , and 
1 /N 3 . The other higher eigenvalues are randomly chosen 
from Ai to 1. We numerically find n c to get an accu- 
racy of Pn c (Xo) > 0.99 (i.e., e < 0.01). Each data is 
averaged over 1000 trials. The data are very well fit- 
ted to In n c = A\nN + B with A ~ 1.130, 2.017, and 
3.002 in each case, which is consistent with our analysis 
of Eq. (f2"12j) . Thus, our algorithm can also be used to esti- 
mate the energy gap between the ground and first-excited 
state (in our case D) [13,|4l[. In Fig.[2jb), we present n c 
with respect to e on a log-log scale. The tests are done for 
three cases: N = 10, 50, and 100. D is fixed to be 1/N. 
Each data are also averaged over 1000 trials. The data 
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FIG. 2: (a) The N versus n c graph is given on a log-log scale. 
Here we consider three cases: D = 1/N, 1/N 2 , and 1/N 3 , 
where D = |Ai — Ao |- Each data is made by averaging 1000 
trials. They are very well fitted to log n c — A log N + B with 
A ~ 1.130 (solid line), 2.017 (dashed line), and 3.002 (dotted 
line) in each case of D. (b) We also present n c versus e graph 
on a log-log scale. We consider three cases: N = 10, 50, and 
100. D is fixed to 1/N. The data are also averaged over 1000 
trials. The data are also well fitted to log n c = A log e+B with 
A ~ -0.187 (solid line), -0.191 (dashed line), and -0.192 
(dotted line) in each case of N. 



are well fitted to logn c = Aloge + B with A ~ -0.187, 
—0.191, and —0.192 for each case. This implies that, for 
the given model Hamiltonian, n c is consistently given in 
the order of C(l/e Q ), where a ~ 0.19 < 1. 

In summary, we have proposed a recursive quantum al- 



gorithm to find the lowest eigenstate of a general Hamil- 
tonian. In our algorithm, a finite number of transfor- 
mations are constructed recursively and applied to an 
arbitrary initial state. Then, one can produce a state ar- 
bitrary close to the lowest eigenstate. Remarkably, our 
algorithm does not depend on any structure of the Hamil- 
tonian. Further, our algorithm does not get trapped in 
the local-minima of the Hamiltonian. The number of re- 
cursions required for the desired accuracy is bounded by 
the order of the difference between the two lowest eigen- 
values. We have numerically investigated the effective- 
ness and efficiency of our algorithm for a simple (but 
non-trivial) hamiltonian model. The results are good 
agreement with our theoretical analyses. We expect that 
our algorithm is widely used for finding lowest eigenstates 
in various specific problems regarding high dimensional 
complex systems. 
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